#include "cfouriertransform.h"
#include <QFile>
#include <QString>
#include <QTextStream>
#include <QRegExp>
#include <cmath>
#include <QDebug>
#include <math.h>
#include <QImage>
#include <QDebug>

QImage *CFourierTransform::GetImg() const
{
    return m_ImgData;
}

QImage *CFourierTransform::GetData() const
{
    return m_ImgData;
}

QImage *CFourierTransform::GetDataData()
{
    return m_Data;
}

void CFourierTransform::SetData(QImage *_data)
{
    m_Data = _data;
}
bool CFourierTransform::IsGrey(){
    return m_Data->allGray();
}

QImage *CFourierTransform::FFT2ImgRGB(Complex** _inputR,Complex** _inputG,Complex** _inputB, bool _outMagnitude)
{
    if (!m_Data)
        return NULL;

    int i, j;
    int nx = m_Width;           // Check this!
    int ny = m_Height;          // Check this!
    float maxR;
     float maxG;
      float maxB;
    float **logFFTR = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        logFFTR[ix] = new float[ny];

    float **logPhaseFFTR = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        logPhaseFFTR[ix] = new float[ny];

    float **FourierValuesR = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        FourierValuesR[ix] = new float[ny];

    float **FourierPhaseR = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        FourierPhaseR[ix] = new float[ny];

    int **normalizedFFTR = new int*[nx];

    for (int ix = 0; ix < nx; ix++)
        normalizedFFTR[ix] = new int[ny];

    int **normalizedPhaseFFTR = new int*[nx];

    for (int ix = 0; ix < nx; ix++)
        normalizedPhaseFFTR[ix] = new int[ny];



    float **logFFTG = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        logFFTG[ix] = new float[ny];

    float **logPhaseFFTG = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        logPhaseFFTG[ix] = new float[ny];

    float **FourierValuesG = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        FourierValuesG[ix] = new float[ny];

    float **FourierPhaseG = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        FourierPhaseG[ix] = new float[ny];

    int **normalizedFFTG = new int*[nx];

    for (int ix = 0; ix < nx; ix++)
        normalizedFFTG[ix] = new int[ny];

    int **normalizedPhaseFFTG = new int*[nx];

    for (int ix = 0; ix < nx; ix++)
        normalizedPhaseFFTG[ix] = new int[ny];





    float **logFFTB = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        logFFTB[ix] = new float[ny];

    float **logPhaseFFTB = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        logPhaseFFTB[ix] = new float[ny];

    float **FourierValuesB = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        FourierValuesB[ix] = new float[ny];

    float **FourierPhaseB = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        FourierPhaseB[ix] = new float[ny];

    int **normalizedFFTB = new int*[nx];

    for (int ix = 0; ix < nx; ix++)
        normalizedFFTB[ix] = new int[ny];

    int **normalizedPhaseFFTB = new int*[nx];

    for (int ix = 0; ix < nx; ix++)
        normalizedPhaseFFTB[ix] = new int[ny];





    uint col;

    if (m_ImgData)
        delete m_ImgData;

    m_ImgData = new QImage(m_Width, m_Height, m_Data->format());

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            FourierValuesR[i][j] = _inputR[i][j].GetMagnitude();
            FourierPhaseR[i][j] = _inputR[i][j].GetPhase();
            logFFTR[i][j] = (float) log(1.0 + 1000.0*FourierValuesR[i][j]);
            logPhaseFFTR[i][j] = (float) log(1 + abs(FourierPhaseR[i][j]));

            FourierValuesG[i][j] = _inputG[i][j].GetMagnitude();
            FourierPhaseG[i][j] = _inputG[i][j].GetPhase();
            logFFTG[i][j] = (float) log(1.0 + 1000.0*FourierValuesG[i][j]);
            logPhaseFFTG[i][j] = (float) log(1 + abs(FourierPhaseG[i][j]));

            FourierValuesB[i][j] = _inputB[i][j].GetMagnitude();
            FourierPhaseB[i][j] = _inputB[i][j].GetPhase();
            logFFTB[i][j] = (float) log(1.0 + 1000.0*FourierValuesB[i][j]);
            logPhaseFFTB[i][j] = (float) log(1 + abs(FourierPhaseB[i][j]));
        }

    m_FourierPhaseR = FourierPhaseR;
    m_FourierPhaseG = FourierPhaseG;
    m_FourierPhaseG = FourierPhaseB;

    // Magnitude
    maxR = logFFTR[0][0];

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            if (logFFTR[i][j] > maxR)
                maxR = logFFTR[i][j];
        }

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            logFFTR[i][j] = logFFTR[i][j] / maxR;
        }

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            normalizedFFTR[i][j] = (int)(255 * logFFTR[i][j]);
        }

    // Magnitude
    maxG = logFFTG[0][0];

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            if (logFFTG[i][j] > maxG)
                maxG = logFFTG[i][j];
        }

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            logFFTG[i][j] = logFFTG[i][j] / maxG;
        }

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            normalizedFFTG[i][j] = (int)(255 * logFFTG[i][j]);
        }

    // Magnitude
    maxB = logFFTB[0][0];

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            if (logFFTB[i][j] > maxB)
                maxB = logFFTB[i][j];
        }

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            logFFTB[i][j] = logFFTB[i][j] / maxB;
        }

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            normalizedFFTB[i][j] = (int)(255 * logFFTB[i][j]);
        }
    // Phase
    logPhaseFFTR[0][0] = 0;
    maxR = logPhaseFFTR[1][1];
    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            if (logPhaseFFTR[i][j] > maxR)
                maxR = logPhaseFFTR[i][j];
        }
    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            logPhaseFFTR[i][j] = logPhaseFFTR[i][j] / maxR;
        }
    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            normalizedPhaseFFTR[i][j] = (int)(255 * logPhaseFFTR[i][j]);
        }

    // Phase
    logPhaseFFTG[0][0] = 0;
    maxG = logPhaseFFTG[1][1];
    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            if (logPhaseFFTG[i][j] > maxG)
                maxG = logPhaseFFTG[i][j];
        }
    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            logPhaseFFTG[i][j] = logPhaseFFTG[i][j] / maxG;
        }
    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            normalizedPhaseFFTG[i][j] = (int)(255 * logPhaseFFTG[i][j]);
        }


    // Phase
    logPhaseFFTB[0][0] = 0;
    maxB = logPhaseFFTB[1][1];
    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            if (logPhaseFFTB[i][j] > maxB)
                maxB = logPhaseFFTB[i][j];
        }
    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            logPhaseFFTB[i][j] = logPhaseFFTB[i][j] / maxB;
        }
    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            normalizedPhaseFFTB[i][j] = (int)(255 * logPhaseFFTB[i][j]);
        }






    // Create img from normal FFT (spectrum from magnitude or spectrum from phase)
    if (_outMagnitude)
    {
        for (i = 0; i <= m_Width - 1; i++)
            for (j = 0; j <= m_Height - 1; j++)
            {
                col = Val2ColRGB(normalizedFFTR[i][j],normalizedFFTG[i][j],normalizedFFTB[i][j]);
                m_ImgData->setPixel(i, j, col);
            }
    }
    else
    {
        for (i = 0; i <= m_Width - 1; i++)
            for (j = 0; j <= m_Height - 1; j++)
            {
                col = Val2ColRGB(normalizedPhaseFFTR[i][j],normalizedPhaseFFTG[i][j],normalizedPhaseFFTB[i][j]);
                m_ImgData->setPixel(i, j, col);
            }
    }

    return m_ImgData;
}


QImage *CFourierTransform::FFT2Img(Complex** _input, bool _outMagnitude)
{
    if (!m_Data)
        return NULL;

    int i, j;
    int nx = m_Width;           // Check this!
    int ny = m_Height;          // Check this!
    float max;
    float **logFFT = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        logFFT[ix] = new float[ny];

    float **logPhaseFFT = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        logPhaseFFT[ix] = new float[ny];

    float **FourierValues = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        FourierValues[ix] = new float[ny];

    float **FourierPhase = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        FourierPhase[ix] = new float[ny];

    int **normalizedFFT = new int*[nx];

    for (int ix = 0; ix < nx; ix++)
        normalizedFFT[ix] = new int[ny];

    int **normalizedPhaseFFT = new int*[nx];

    for (int ix = 0; ix < nx; ix++)
        normalizedPhaseFFT[ix] = new int[ny];

    uint col;

    if (m_ImgData)
        delete m_ImgData;

    m_ImgData = new QImage(m_Width, m_Height, m_Data->format());

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            FourierValues[i][j] = _input[i][j].GetMagnitude();
            FourierPhase[i][j] = _input[i][j].GetPhase();
            logFFT[i][j] = (float) log(1.0 + 1000.0*FourierValues[i][j]);
            logPhaseFFT[i][j] = (float) log(1 + abs(FourierPhase[i][j]));
        }

    m_FourierPhase = FourierPhase;

    // Magnitude
    max = logFFT[0][0];

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            if (logFFT[i][j] > max)
                max = logFFT[i][j];
        }

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            logFFT[i][j] = logFFT[i][j] / max;
        }

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            normalizedFFT[i][j] = (int)(255 * logFFT[i][j]);
        }
    // Phase
    logPhaseFFT[0][0] = 0;
    max = logPhaseFFT[1][1];
    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            if (logPhaseFFT[i][j] > max)
                max = logPhaseFFT[i][j];
        }
    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            logPhaseFFT[i][j] = logPhaseFFT[i][j] / max;
        }
    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            normalizedPhaseFFT[i][j] = (int)(255 * logPhaseFFT[i][j]);
        }

    // Create img from normal FFT (spectrum from magnitude or spectrum from phase)
    if (_outMagnitude)
    {
        for (i = 0; i <= m_Width - 1; i++)
            for (j = 0; j <= m_Height - 1; j++)
            {
                col = Val2Col(normalizedFFT[i][j]);
                m_ImgData->setPixel(i, j, col);
            }
    }
    else
    {
        for (i = 0; i <= m_Width - 1; i++)
            for (j = 0; j <= m_Height - 1; j++)
            {
                col = Val2Col(normalizedPhaseFFT[i][j]);
                m_ImgData->setPixel(i, j, col);
            }
    }

    return m_ImgData;
}



uint CFourierTransform::Val2Col(int _val)
{
    uint retCol;
    int red;
    int green;
    int blue;

    blue = _val & 0xFF;
    red = blue;
    green = red;

    if (red > 255)
        red = 255;
    else if (red < 0)
        red = 0;

    if (green > 255)
        green = 255;
    else if (green < 0)
        green = 0;

    if (blue > 255)
        blue = 255;
    else if (blue < 0)
        blue = 0;

    retCol = 0xFF;
    retCol <<= 8;
    retCol |= red;
    retCol <<= 8;
    retCol |= green;
    retCol <<= 8;
    retCol |= blue;

    return retCol;
}
uint CFourierTransform::Val2ColRGB(int _valR, int _valG, int _valB)
{
    uint retCol;
    int red;
    int green;
    int blue;

    blue = _valB & 0xFF;
    red = _valR & 0xFF;
    green = _valG & 0xFF;

    if (red > 255)
        red = 255;
    else if (red < 0)
        red = 0;

    if (green > 255)
        green = 255;
    else if (green < 0)
        green = 0;

    if (blue > 255)
        blue = 255;
    else if (blue < 0)
        blue = 0;

    retCol = 0xFF;
    retCol <<= 8;
    retCol |= red;
    retCol <<= 8;
    retCol |= green;
    retCol <<= 8;
    retCol |= blue;

    return retCol;
}

void CFourierTransform::FillTestFFTData()
{
    QFile file("output.txt");
    file.open(QIODevice::ReadOnly | QIODevice::Text);

    if (!file.isOpen())
        return;

    QTextStream in(&file);
    QString line;
    QRegExp rxlen("^(.*):(.*)$");
    int pos;
    double _re;
    double _im;
    QString _imStr;
    QString _reStr;

    m_TestFFTData = new Complex*[m_Width]();


    for (int i = 0; i < m_Width; i++)
    {
        m_TestFFTData[i] = new Complex[m_Height];

        for (int j = 0; j < m_Height; j++)
        {
            line = in.readLine();
            pos = rxlen.indexIn(line);

            if (pos > -1)
            {
                _imStr = rxlen.cap(1);
                _reStr = rxlen.cap(2);
                _im = rxlen.cap(1).toDouble();
                _re = rxlen.cap(2).toDouble();
                m_TestFFTData[i][j].im = _im;
                m_TestFFTData[i][j].re = _re;
            }
        }
    }
}

QImage* CFourierTransform::PhaseMod2Img(float** _phase)
{
    int i, j;
    int nx = m_Width;
    int ny = m_Height;
    float max;
    float **logPhaseFFT = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        logPhaseFFT[ix] = new float[ny];

    float **FourierPhase = new float*[nx];

    for (int ix = 0; ix < nx; ix++)
        FourierPhase[ix] = new float[ny];

    int **normalizedPhaseFFT = new int*[nx];
    for (int ix = 0; ix < nx; ix++)
        normalizedPhaseFFT[ix] = new int[ny];

    uint col;

    if (m_ImgData)
        delete m_ImgData;

    m_ImgData = new QImage(m_Width, m_Height, m_Data->format());

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            logPhaseFFT[i][j] = (float) log(1 + abs(_phase[i][j]));
        }

    // Phase
    logPhaseFFT[0][0] = 0;
    max = logPhaseFFT[1][1];

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            if (logPhaseFFT[i][j] > max)
                max = logPhaseFFT[i][j];
        }
    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            logPhaseFFT[i][j] = logPhaseFFT[i][j] / max;
        }
    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            normalizedPhaseFFT[i][j] = (int)(255 * logPhaseFFT[i][j]);
        }

    for (i = 0; i <= m_Width - 1; i++)
        for (j = 0; j <= m_Height - 1; j++)
        {
            col = Val2Col(normalizedPhaseFFT[i][j]);
            m_ImgData->setPixel(i, j, col);
        }

    return m_ImgData;
}
void CFourierTransform::ImgToComplex()
{
    uint col;
    int red;
    int green;
    int blue;
    int sumCols;

    fft = new Complex*[m_Width]();
    for (int i = 0; i < m_Width; i++)
    {
        fft[i] = new Complex[m_Height];

        for (int j = 0; j < m_Height; j++)
        {
            col = m_Data->pixel(i, j);
            red = (col >> 16) & 0xFF;
            green = (col >> 8) & 0xFF;
            blue = col & 0xFF;
            sumCols = (red + green + blue) / 3;
            col = 0xFF;
            col <<= 8;
            col |= sumCols;
            col <<= 8;
            col |= sumCols;
            col <<= 8;
            col |= sumCols;
            fft[i][j].re = (double) sumCols;
            fft[i][j].im = 0;
        }
    }
}
void CFourierTransform::ImgOutToComplex()
{
    uint col;
    int red;
    int green;
    int blue;
    int sumCols;

    forIFFT = new Complex*[m_Width]();
    for (int i = 0; i < m_Width; i++)
    {
        forIFFT[i] = new Complex[m_Height];

        for (int j = 0; j < m_Height; j++)
        {
            col = m_ImgData->pixel(i, j);
            red = (col >> 16) & 0xFF;
            green = (col >> 8) & 0xFF;
            blue = col & 0xFF;
            sumCols = (red + green + blue) / 3;
            col = 0xFF;
            col <<= 8;
            col |= sumCols;
            col <<= 8;
            col |= sumCols;
            col <<= 8;
            col |= sumCols;
            forIFFT[i][j].re = (double) sumCols;
            forIFFT[i][j].im = 0;
        }
    }
}
void CFourierTransform::ImgToComplexRGB()
{
    uint col;
    int red;
    int green;
    int blue;
    int sumCols;

    fftR = new Complex*[m_Width]();
    fftG = new Complex*[m_Width]();
    fftB = new Complex*[m_Width]();
    for (int i = 0; i < m_Width; i++)
    {
        fftR[i] = new Complex[m_Height];
        fftG[i] = new Complex[m_Height];
        fftB[i] = new Complex[m_Height];

        for (int j = 0; j < m_Height; j++)
        {
            col = m_Data->pixel(i, j);
            red = (col >> 16) & 0xFF;
            green = (col >> 8) & 0xFF;
            blue = col & 0xFF;
            //sumCols = (red + green + blue) / 3;
            col = 0xFF;
            col <<= 8;
            col |= red;
            col <<= 8;
            col |= green;
            col <<= 8;
            col |= blue;
            fftR[i][j].re = (double) red;
            fftR[i][j].im = 0;
            fftG[i][j].re = (double) green;
            fftG[i][j].im = 0;
            fftB[i][j].re = (double) blue;
            fftB[i][j].im = 0;
        }
    }
}
Complex** CFourierTransform::ShiftFFT(Complex **tab){
    int i,j;
    Complex**shiftedFFT = new Complex*[m_Width]();
    for(i=0; i<m_Width;i++)
        shiftedFFT[i] = new Complex[m_Height];
    for(i=0; i<=(m_Width/2)-1;i++){
        for(j=0;j<=(m_Height/2)-1; j++){
            shiftedFFT[i + (m_Width / 2)][j + (m_Height/2)] = tab[i][j];
            shiftedFFT[i][j] = tab[i + (m_Width / 2)][j + (m_Height / 2)];
            shiftedFFT[i + (m_Width/2)][j] = tab[i][j + (m_Height / 2)];
            shiftedFFT[i][j + (m_Width / 2)] = tab[i + (m_Width / 2)][j];
        }
    }
    return shiftedFFT;
}

void CFourierTransform::IFFT(){
    m_Width = m_Data->width();
    m_Height = m_Data->height();
    uint col;
    if(m_Data->allGray()){
        //ImgToComplex();
        Complex **tmp = ShiftFFT(fft);
        FFT2D(tmp, -1);
        if (m_ImgData)
            delete m_ImgData;

        m_ImgData = new QImage(m_Width, m_Height, m_Data->format());
        for (int i = 0; i < m_Width ; i++)
        {
            for (int j = 0; j < m_Height ; j++)
            {
                col = Val2Col((int)( tmp[i][j].GetMagnitude()));
                m_ImgData->setPixel(i, j, col);
            }
        }
    }
    else{
        //ImgToComplexRGB();
        Complex **tmpR = ShiftFFT(fftR);
        Complex **tmpG = ShiftFFT(fftG);
        Complex **tmpB = ShiftFFT(fftB);
        FFT2D(tmpR, -1);
        FFT2D(tmpG, -1);
        FFT2D(tmpB, -1);
        if (m_ImgData)
            delete m_ImgData;

        m_ImgData = new QImage(m_Width, m_Height, m_Data->format());
        for (int i = 0; i < m_Width ; i++)
        {
            for (int j = 0; j < m_Height ; j++)
            {
                col = Val2ColRGB((int)tmpR[i][j].GetMagnitude(),(int)tmpG[i][j].GetMagnitude(),(int)tmpB[i][j].GetMagnitude());
                m_ImgData->setPixel(i, j, col);
            }
        }
    }
}
void CFourierTransform::FFT(){
    m_Width = m_Data->width();
    m_Height = m_Data->height();
    if(m_Data->allGray()){
        qDebug()<<"szary";
        ImgToComplex();
        FFT2D(fft, 1);
        fft=ShiftFFT(fft);
    }
    else{
        qDebug()<<"kolorowy";
        ImgToComplexRGB();
        FFT2D(fftR, 1);//nie trzeba podawać innych rzeczy niiż kierunek bo inne sa globalne
        FFT2D(fftG, 1);
        FFT2D(fftB, 1);
        fftR = ShiftFFT(fftR);
        fftG = ShiftFFT(fftG);
        fftB = ShiftFFT(fftB);
    }
}

void CFourierTransform::FFT2D(Complex** complex, int dir){

    Complex *buf = new Complex[m_Width];
    int i,j;
    int mm;
            //nie wiem po co to jest i daje z dupy wynki
    //fft[(int)(m_Width/2)][(int)(m_Height/2)].re=0.0;
    //fft[(int)(m_Width/2)][(int)(m_Height/2)].im=0.0;
    /* Transform the rows*/
    for (j=0;j<m_Height;j++) {
      for (i=0;i<m_Width;i++) {
          buf[i].re = complex[i][j].re;
          buf[i].im = complex[i][j].im;

      }
      mm = log(m_Width)/log(2);
      FFT1D(dir, m_Width, mm, buf);

      for (i=0;i<m_Width;i++) {
          complex[i][j].re = buf[i].re;
          complex[i][j].im = buf[i].im;
      }
    }

    /* Transform the columns*/
       delete buf;
       Complex *buf1 = new Complex[m_Height];

    for (i=0;i<m_Width;i++) {
      for (j=0;j<m_Height;j++) {
          buf1[j].im = complex[i][j].im;
          buf1[j].re = complex[i][j].re;
      }
      mm = log(m_Height)/log(2);
      FFT1D(dir,m_Height, mm, buf1);
      for (j=0;j<m_Height;j++) {
          complex[i][j].im = buf1[j].im;
          complex[i][j].re = buf1[j].re;
      }
    }


    return;
}


void CFourierTransform::FFT1D(int dir, int m, int mm, Complex *buf){
    double tmpRe, tmpIm;
    long i,i1,i2,j,k,l,l1,l2;
    double u1,u2,z;
    /* reverse-binary reindexing */
    //m-ilość punktów
    i2 = m >> 1;
    j = 0;

    for (i = 0; i < m-1; i++) {
      if (i < j) {
         tmpRe = buf[i].re;
         tmpIm = buf[i].im;
         buf[i].re = buf[j].re;
         buf[i].im = buf[j].im;
         buf[j].re = tmpRe;
         buf[j].im = tmpIm;
      }
      k = i2;
      while (k <=j) {
         j -= k;
         k >>= 1;
      }
      j += k;
    }

    /* Compute the FFT */

    double c1 = -1.0;
    double c2 = 0.0;
    l2=1;
    for(l=0;l<mm;l++){
        l1 = l2;
        l2 <<= 1;
        u1 = 1.0;
        u2 = 0.0;

        for(j=0; j<l1; j++){
           for(i=j; i<m; i+=l2){
               i1=i+l1;
               tmpRe = u1 * buf[i1].re - u2*buf[i1].im;
               tmpIm = u1 * buf[i1].im + u2*buf[i1].re;

               buf[i1].re = buf[i].re - tmpRe;
               buf[i1].im = buf[i].im - tmpIm;
               buf[i].re += tmpRe;
               buf[i].im += tmpIm;
           }
           z = u1 * c1 - u2 * c2;
           u2 = u1*c2 + u2 *c1;
           u1=z;
        }
    c2 = sqrt((1.0 - c1)/2.0);
    c1 = sqrt((1.0 + c1)/2.0);

    if(dir == 1)
       c2 = -c2;
    }
   /* Scaling for forward transform */
   if (dir == 1) {
      for (i=0;i<m;i++) {
         buf[i].im /= (double)m;
         buf[i].re /= (double)m;
      }
   }

   return;
}
/*void CFourierTransform::FFT2DNew(Complex** complex, int dir){

    Complex *buf = new Complex[m_Width];
    int i,j;

    SymetricFFT();


    for (j=0;j<m_Height;j++)
    {
        for (i=0;i<m_Width;i++)
        {
            buf[i].re = complex[i][j].re;
            buf[i].im = complex[i][j].im;

        }
    }

    FFT1DNew(buf, buf, m_Width, dir);

    for (j=0;j<m_Height;j++)
        for (i=0;i<m_Width;i++)
        {
            complex[i][j].re = buf[i].re;
            complex[i][j].im = buf[i].im;
        }

    delete buf;
    Complex *buf1 = new Complex[m_Height];

    for (i=0;i<m_Width;i++)
    {
        for (j=0;j<m_Height;j++)
        {
            buf1[j].im = complex[i][j].im;
            buf1[j].re = complex[i][j].re;
        }
    }

    FFT1DNew(buf1, buf1, m_Height, dir);

    for (i=0;i<m_Width;i++)
        for (j=0;j<m_Height;j++)
        {
            complex[i][j].im = buf1[j].im;
            complex[i][j].re = buf1[j].re;
        }

    return;
}
*/
/*
void CFourierTransform::SymetricFFT()
{
    m_iN = m_Width;
    m_ppW = new Complex*[m_iN + 1];

    for (int i = 0; i <= m_iN; i++)
    {
        m_ppW[i] = new Complex[m_iN + 1];

        for (int j = 0; i > 0 && j <= m_iN; j++)
        {
            float fPs = (m_f2Pi/i) * j;
            m_ppW[i][j].re = cos(fPs);
            m_ppW[i][j].im = -sin(fPs);
        }
    }
}

void CFourierTransform::FFT1DNew(Complex* F, Complex* f, int iN, int iFact)
{
    int iM = iN / 2;

    if (iM == 1)
    {
        F[0] = (f[0] + f[1]).DivByNumber(m_fSqrt2);
        F[1] = (f[0] - f[1]).DivByNumber(m_fSqrt2);
    }
    else
    {
        Complex* g = new Complex[iM];
        Complex* h = new Complex[iM];
        Complex* G = new Complex[iM];
        Complex* H = new Complex[iM];

    for (int x = 0; x < iM; x++)
    {
        g[x] = f[2*x];
        h[x] = f[2*x+1];
    }

        FFT1DNew(G, g, iM, iFact);
        FFT1DNew(H, h, iM, iFact);

    for (int u = 0; u < iM; u++)
    {
        Complex xW(m_ppW[iN][u].GetMagnitude(), iFact * m_ppW[iN][u].GetPhase());
        Complex xAux = H[u] * xW;
        F[u]   = (G[u] + xAux).DivByNumber(m_fSqrt2);
        F[iM+u] = (G[u] - xAux).DivByNumber(m_fSqrt2);
    }

        delete [] g;
        delete [] h;
        delete [] G;
        delete [] H;

    }
}
*/
